####################### Reproduce Appendix Figure 2 (SRBI)###############################

### PREP
  rm(list = ls())
  #setwd(..)
  smoking <- read.dta("srbixx.dta")
  library(gam)
  library(Hmisc)

### RELATIVE VS. ATTRIBUTABLE RISK 
  # Quadrant 2 graph - page 2 of pdf output 
  # Quadrant 3 graph - page 3 of pdf output 

  smoking.gam1 <- gam(currsmok2 ~ s(rd) + s(rr),
                     subset = control==1&rd>-500,
                     data=smoking, family=gaussian(), weights=weight1)
  smoking.gam1 <- gam(currsmok2 ~ s(rd) + s(rr),
                     subset = control==1&rd>0&rd<800,
                     data=smoking, family=gaussian(), weights=weight1)
  pdf("srbi_attributable.pdf")
    plot(smoking.gam1, xlab="Relative Risk", se=TRUE, ylab="Probability of Being a Current Smoker", ask=F, ylim=c(-.4,.2))
    plot(smoking.gam1, xlab="Attributable Risk", se=TRUE, ylab="Probability of Being a Current Smoker", ask=F)
  dev.off()

### RELATIVE VS. ABSOLUTE RISK 
  # Quadrant 1 graph - page 2 of pdf output
  # Quadrant 4 graph - page 3 of pdf output 
  
  smoking.gam1 <- gam(currsmok2 ~ s(pq52) + s(rr),
                     subset = control==1&rd>-500,
                     data=smoking, family=gaussian(),weights=weight1)
  
  smoking.gam1 <- gam(currsmok2 ~ s(pq52) + s(rr),
                     subset = control==1&rd>-500&rd>0&rd<1000,
                     data=smoking, family=gaussian(),weights=weight1)
  pdf("srbi_absolute.pdf")
    plot(smoking.gam1, xlab="Relative Risk", se=TRUE, ylab="Probability of Being a Current Smoker", ask=F)
    plot(smoking.gam1, xlab="Absolute Risk", se=TRUE, ylab="Probability of Being a Current Smoker", ask=F)
  dev.off()

  